######### This script will summarize fragmentation statistics for 7 species of skink

#### Establish base directory

base.dir = '/home1/99/jc152199/MAXENT/'
setwd(base.dir)

#### Load library

library('SDMTools')

##### Need to loop through ASCIIs for multiple species
##### First identify ASCII directory

ascii.dir = '/home1/99/jc152199/MAXENT/output/'

#### Directory to write files to

out.dir = '/home1/99/jc152199/FirstPublication/SkinkFragStats/'

##### Identify species

spp = dir(ascii.dir)

#### Remove directories which aren't skinks

spp = spp[c(1,13:15,17:19)]

###### Now loop through and identify individual ASCII for each species

for (s in spp)

	{
	
	#### Species loop
	
	for (m in c('microCLIM','BIOCLIM'))
	
		{
		
		#### Model loop
		
		#### Read in an ASCII of realized distribution
		
		tasc = read.asc.gz(paste(ascii.dir,s,'/',m,'/output/',s,'_',m,'_Realized.asc.gz',sep=''))
		
		#### Keep an unchanged copy of tasc called base.asc
		
		base.asc =  tasc
		
		#### Convert to a binary matrix
		#### Don't need to apply MaxEnt thresh because that was done before producing the realized distribution
		
		tasc[which(tasc>0)]=1
		
		#### Run class-stat on tasc to determine total area for the species distribution
		
		cstat = ClassStat(tasc,cellsize=250,bkgd=NA,latlon=TRUE)[2,]
		
		#### Calculate summed environmental suitability across the whole distribution
		#### This is a proxy for abundance
		
		cstat$ESsum = sum(base.asc,na.rm=T)
		
		### Remove non-important stats
		
		cstat = cstat[,c(2,3,34,37,39)]
		
		### Write to the species specific directory
		
		write.csv(cstat, file=paste(out.dir,s,'/',m,'_ClassStat.csv',sep=''), row.names=F)
		
		### Now determine stats for individual patches
		#### Run ConnCompLabel to identify individual patches
		
		pasc = ConnCompLabel(tasc)
		
		#### Run Patch Stat on pasc
		
		pstat = PatchStat(pasc, cellsize=250, latlon=TRUE)
		
		### Subset pstat to only stats of interest
		### Removing the first row (which describes the shape of the un-occupied patch)
		
		pstat = pstat[-1,c(1,6,9)]
		
		### Now need to calculate summed ES for each patch
		### Use a for loop to cycle through individual patches identified in pstat
		
		### Create a blank column in pstat first to bind data to
		
		pstat$ESsum = NA
		
		### Position tracker
		
		i=1
		
		for (p in pstat$patchID)
		
			{
			
			### Identify row/col positions in pasc that match p
			
			ppos = as.data.frame(which(pasc==p, arr.ind = T))
			
			### Identify the corresponding geographic positions
			
			ppos$lat = getXYcoords(pasc)$x[ppos$row]
			ppos$long = getXYcoords(pasc)$y[ppos$col]
			
			### Extract data from base.asc (which contains ES values) at the above positions
			
			pES = extract.data(cbind(ppos$lat,ppos$long),base.asc)
		
			### Sum the values in pES to get the total environmental suitability (potential abundance) for that patch
			
			pstat$ESsum[which(pstat$patchID==p)] = sum(pES)
			
			### Report progress
			
			cat('\n',(i/length(pstat$patchID))*100,' - Percent Complete',sep='')
			
			### Change position tracker
			
			i=i+1
			
			}
			
		### Close loop
			
		### Write out pstat to the out.dir
		
		write.csv(pstat, file=paste(out.dir,s,'/',m,'_PatchStat.csv',sep=''), row.names=F)
		
		}
		
	cat('\n',s,' - Completed\n',sep='')	
		
	}
	
#### Done
